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Abstract 

In this paper we present a systematic and rigorous method for calculating the 
diffusion tensor for a Brownian particle moving in a periodic potential which is 
valid in arbitrary dimensions and for all values of the dissipation. We use this 
method to obtain an explicit formula for the diffusion coefficient in one dimension 
which is valid in the underdamped limit, and we also obtain higher order correc- 
tions to the Lifson- Jackson formula for the diffusion coefficient in the overdamped 
limit. A numerical method for calculating the diffusion coefficient is also devel- 
oped and is shown to perform extremely well for all values of the dissipation. 



1 Introduction 

Brownian motion in periodic and random potentials has been a very active area of re- 
search for many decades. Apart from the well established applications to electronics 
ll34l l26l and to solid state physics such as superionic conductors, the Josephson tun- 
neling junction (T) and surface diffusion [ 10|, new and exciting applications to physics 
(self-assembled molecular film growth, catalysis, surface-bound nanostructures) and to 
biology (stochastic modeling of molecular and Brownian motors [27 1) keep the subject 
of Brownian motion at the forefront of current research, both theoretical and experi- 
mental. 

Despite the fact that Brownian motion in periodic potentials has been studied exten- 
sively f30l Ch. 11], and many analytical and numerical results have been obtained, 
there are still many open questions, in particular in the underdamped, multidimensional 
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case. The main purpose of this paper is to develop a general method for calculating the 
diffusion tensor D of the Brownian particle in arbitrary dimensions, and to then use 
this method for setting up an efficient numerical method for computing D. Further- 
more, we will show that our method for calculating D will enable us to study various 
asymptotic limits of physical interest in a systematic and rigorous fashion. 

The dynamics of a Brownian particle moving in a periodic potential is governed by 
the Langevin equation 

x = -W(x) - 7± + £, (1.1) 

where x(t) denotes the particle position, V(x) is a smooth periodic potential, 7 denotes 
the friction coefficient and £ (t) is a white noise Gaussian process with correlation func- 
tion 

&(*)&(*)) = 2lk B T5 tJ 6(t ~ s), i,j = l,...,d. 

where kg is Boltzmann's constant and T is the absolute temperature, in accordance 
with the fluctuation-dissipation theorem. We use the notation (•) to denote ensemble 
average. We will also write /3 = (fcsT)" 1 and = y^y^W, where W(t) is 
a standard Brownian motion in M. d . Throughout this work we will assume that the 
diffusing particle is of unit mass, rn = 1. 

The Langevin equation ( 11.11 ) has been studied extensively as a theoretical model for 
the diffusion of adsorbates on crystal surfaces |[T0ll2TI . In this setting, q(t) represents 
the position of the diffusing particle, V(q) the substrate potential and the friction and 
noise terms represent the interaction of the diffusing particle with the phonon heat 
bath ED 

It is well known that at low temperatures (which is usually the regime of physical 
interest), the diffusing particle performs a hopping motion (random walk) between the 
local minima of the potential. This hopping motion is characterized by the mean square 
jump length (I 2 ) and the hopping rate k (or, equivalently, the mean escape time r with 
k = i). These two quantities are related to the diffusion coefficient through the 
formula (in one dimension) 

D = (P)k. (1.2) 

Knowledge of two of the three quantities (D, (£ 2 ), n) is sufficient for the calculation 
of the third. Of course, the diffusion tensor can also be defined either in terms of the 
ensemble average of the second moment, i.e. 

D = Km I<( ? (t) - q(0)) ® (q(t) - q(0))), (1.3) 

or in terms of the time integral of the velocity autocorrelation function, i.e. through the 
Green-Kubo formula 

poo 

D= (p(t) <g> p(0)) dt (1.4) 
Jo 

with p(t) = q(t). In the above formulas ® stands for the tensor product between two 
vectors in W l . We remark that, although formulas d 1 -3b and dl.4l ) are valid in arbitrary 
dimensions, it is not clear how to interpret formula ( 11.21 ) in dimensions higher than 1. 

Of particular interest is the dependence of the diffusion coefficient D as well as the 
jump rate k and (mean square) jump length (i 2 ) on the friction coefficient. In particular, 
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it is well known that in the underdamped regime the particle diffusion is dominated by 
the occurrence of long jumps J2] [4] QT). Apart from theoretical investigations and 
numerical simulations based on the Langevin dynamics ( 11.11 ), the occurrence of long 
jumps in the underdamped regime is also verified by means of molecular dynamics 
simulations of a model for CO/Ni(lll) (7) and is also, by now, a well established 
experimental result [32 21]. Indeed, in the case of weak adsorbate-substrate interaction 
i.e. in the case of weak coupling between the diffusing particle and the heat bath 
which corresponds to the underdamped dynamics regime, the diffusion mechanism is 
controlled by long jumps, spanning multiple lattice spacings [ 32 , 2 1 ] . The rigorous and 
systematic study of the diffusion process in the weak dissipation regime is still a major 
challenge for theoreticians, in particular in dimensions higher than 1. 

The problem of diffusion in periodic potentials is well studied in one dimension ll30l 
12 [2T) or for separable potentials in two and three dimensions lfT~8l . Kramers' the- 
ory (17| applied to Brownian motion in a periodic potential or the mean first passage 
time method [33 1 enables us to calculate the rate of escape from a local minimum of 
the potential (hopping rate). In the overdamped limit this is sufficient to calculate the 
diffusion coefficient, since only single jumps occur and consequently (i 2 ) = L 2 where 
L is the period of the potential. This leads to the well known Lifson-Jackson formula 
for the diffusion coefficient HU (25] Ch. 13] 

D= Ml , (i. 5) 

eP v W dq f Q 

where = fc^T and Dq denotes the diffusion coefficient of the free particle 

7 

Kramers' formula for the rate of escape enables us to obtain a formula for the dif- 
fusion coefficient which is valid in the moderate-to-strong friction regime, for small 
temperatures [3 > 1 lfTTl[T6l : 

(1.6) 

whereE'f, = V(q M Ax) - V(q M m), = V"(q M iN), u\ = \V"(q M Ax)\- In the 
overdamped, j/u>b 1, small temperature ^ > 1, limit this formula reduces to the 
small temperature asymptotics of equation ( 11.5b : 

D= <W*& e -M (1.7) 

2-7T7 

The calculation of the diffusion coefficient in the underdamped limit requires two cal- 
culations, that of the hopping rate and that of the mean squared jump length (i 2 ). Such 
a calculation was presented in l3Tl[T8l for the case of a cosine potential. For this po- 
tential, a formula for the diffusion coefficient which is valid in the regime 7 <C 1 was 
also obtained by Risken and presented in his monograph (3D): 

1 IT O o 

D =-,W 5 (L8) 
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In contrast to the one dimensional problem, a similar theory in higher dimensions is still 
lacking, except for the overdamped limit. In this limit approximate analytical results for 
certain two-dimensional potentials have been derived in the literature [3 1. Furthermore, 
it is possible to prove, using rigorous mathematical analysis, that the diffusion tensor 
scales like - when 7> 1 lfl2l . 

On the contrary, it is still not clear what the scaling of the (trace of the) diffu- 
sion tensor with the friction constant is in the underdamped, multidimensional case. 
Numerical experiments [2] suggest that this scaling depends crucially on the detailed 
properties of the periodic potential; however, a rigorous and systematic theory for ex- 
plaining the dependence of the diffusion coefficient on the strength of the dissipation 
in arbitrary dimensions is still lacking. 

The situation becomes even more unclear when an external driving force (either 
constant or periodic in time) is present. Brownian motion in tilted periodic potentials 
has only been studied in one dimension [ 30 1 and explicit formulas are only valid in 
the overdamped limit l28l l20l l24l . Furthermore, numerical experiments ll3~5l seem to 
indicate that stochastic resonance in periodic potentials in only possible in dimensions 
higher than one. 

In view of the ubiquity of diffusive motion in periodic potentials in applications, 
it seems to be important to develop the multidimensional theory in a systematic and 
rigorous way, taking into account external driving forces. This paper is a contribution 
towards this goal, and it is a part of our research program on the study of Brownian 
motion in periodic and random potentials 11231 [T2l l24l . 

All studies on the problem of Brownian motion in a periodic potential that have 
been reported in the physics literature rely heavily on the analysis of the Fokker- 
Planck (Kramers-Chandrashekhar) equation which governs the evolution of the transi- 
tion probability density for the Brownian particle. For example, the continued fraction 
expansion method is used in order to solve the Fokker-Planck equation j5] [30) in a 
semi-analytic fashion and to calculate quantities such as the mobility, the intermediate 
structure function and the dynamic structure function [6 1. Or, Kramers' theory is being 
used, which again relies on the study of the Fokker-Planck equation. 

On the other hand, many tools from stochastic analysis [29], the theory of limit 
theorems for Markov processes [8 1 and the emerging field of multiscale analysis ll25l 
are appropriate for the study of this problem and, yet, they have received very little 
attention in the physics community. The purpose of this paper is to use multiscale 
methods such as homogenization theory and singular perturbation theory in order to 
offer a new insight into the problem of Brownian motion in periodic potentials. In par- 
ticular, we derive in a rigorous and systematic way a formula for the diffusion tensor 
D of a Brownian particle moving in a periodic potential in arbitrary dimensions and 
then we use in order to develop an efficient numerical method for calculating D. This 
numerical method is related to the continued fraction expansion method, but is easier to 
implement and to analyze. As a byproduct of our analysis, we derive rigorously a for- 
mula for the diffusion coefficient which is valid in the weak friction limit; furthermore, 
we also calculate higher order corrections to the large 7 asymptotics of the diffusion 
coefficient. 

The rest of the paper is organized as follows. In Section [2] we present the multi- 
scale analysis and we derive a formula for the diffusion coefficient. We also show the 
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equivalence between our formula and the Green-Kubo formula ( 11.4b . In Section (0 we 
derive formulas which are valid in the 7^0 and 7 — * +00 limits. In Section [4] we 
develop the numerical method and we compare the numerical results obtained using 
our method with results obtained from Monte Carlo simulations, from approximate an- 
alytical formulas and from the numerical implementation of formula J1.2t . Section[5] 
is reserved for conclusions. 

2 Multiscale Analysis 

In this section we use multiscale analysis l25l to derive a formula for the diffusion 
tensor of a Brownian particle moving in a periodic potential in arbitrary dimensions. 
We then show the equivalence between this formula and the Green-Kubo formula for 
the diffusion tensor. 

2.1 Derivation of Formula for the Diffusion Tensor 

We start by rescaling the Langevin equation fll.lt 

x = F(x) - 7± + y/2 1 (3- 1 W, (2.1) 

where we have set F(x) — — W(x). We will assume that the potential is periodic 
with period 2tt in every direction. Since we expect that at sufficiently long length and 
time scales the particle performs a purely diffusive motion, we perform a diffusive 
rescaling to the equations of motion ( II . lb : t — > t/e 2 , x — > —. Using the fact that 
W{ct) = -j=W(t) in law we obtain: 

e 2 x = -F (-) - >yx + y/2<y(3- 1 W, 
Introducing p — ex and q = x/e we write this equation as a first order system: 

x = \p, 

V = ^F(q)~^ 1P +^p- 1 W, (2.2) 

<i = jsp, 

with the understanding that q 6 [— it, ir] d and x, p £ R d . Our goal now is to eliminate 
the fast variables p, q and to obtain an equation for the slow variable x. We shall 
accomplish this by studying the corresponding backward Kolmogorov equation using 
singular perturbation theory for partial differential equations. 
Let 

u e (p,q,x,t) =Ef(p(t),q(t),x(t)\p(0) =p,q(0) = q,x(0) = x), 

where E denotes the expectation with respect to the Brownian motion Wif) in the 
Langevin equation and / is a smooth function^ The evolution of the function u e (p, q,x,t) 

1 In other words, we have that 

u e (p,q,x,t) = / f(x,v,t;p,q)p(x,v,t;p,q)[i(p,q)dpdqdxdv, 
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is governed by the backward Kolmogorov equation associated to equations ( 12.21 ) is [25 

= -V ■ ^ x u e + ^ ( - V g V(q) • V p + p • V 9 + 7 ( - p • V p + /8 _1 Ap)) 



"a* 



^£ + ^i)" £ , (2.3) 



where: 



Co = -V q V(q)-V p +p-V q +j(- P -V p + f3- 1 A p ), 
A = p-V x 

The invariant distribution of the fast process {q(t), p(t)} in {—ir,ir] d x M. d is the 
Maxwell-Boltzmann distribution 

p p (q,p) = Z- 1 e-" H < w ') 1 Z = / e -/"T(wO dqdp, 

J[-7T,Tr] d XR t! 

where (g, p) = \ \p\ 2 + V(q). Indeed, we can readily check that 

£oP/3(q,p) = o, 

where Cq denotes the Fokker-Planck operator which is the L 2 -adjoint of the generator 
of the process Co: 

£*of- = V q V(q) ■ Vp/ - p ■ V q f + 7 (Vp • (pf) + /T 1 Ap/) . 

The null space of the generator Co consists of constants in q, p. Moreover, the equation 

- C f = .9, (2.4) 

has a unique (up to constants) solution if and only if 

(g)p:= g{q 1 p)p l3 (q,p)dqdp = 0. (2.5) 



•7T.7T "X.- 



Equation (12.4-b is equipped with periodic boundary conditions with respect to z and is 
such that 

/ \f\ 2 ppdqdp < Co. (2.6) 

J[-Tr,ir] d xR d 

These two conditions are sufficient to ensure existence and uniqueness of solutions (up 
to constants) of equation ( |2~4l > |[T2l [T3l l22l . 



where p(x, v,t;p, q) is the solution of the Fokker-Planck equation and fi(p, q) is the initial distribution. 

2 it is more customary in the physics literature to use the forward Kolmogorov equation, i.e. the Fokker- 
Planck equation. However, for the calculation presented below, it is more convenient to use the backward 
as opposed to the forward Kolmogorov equation. The two formulations are equivalent. See 1231 Ch. 6] for 
details. 
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We assume that the following ansatz for the solution u e holds: 

u £ — it + eui + e 2 u 2 + ... (2.7) 

with Ui = Ui(p,q,x,t), i — 1,2,... being 2tt periodic in q and satisfying condi- 
tion (12.6b . We substitute ( 12.7b into ( 12.3b and equate equal powers in e to obtain the 
following sequence of equations: 

C Q u Q = 0, (2.8a) 
Cqux = -CiUq, (2.8b) 

C u 2 = -A + (2.8c) 
at 

From the first equation in J2.81 l we deduce that u — u (x,t), since the null space of 
Cq consists of functions which are constants in p and q. Now the second equation in 
( 12.81 ) becomes: 

C Q ui = -p ■ V x u Q . 

Since (p) = 0, the right hand side of the above equation is mean-zero with respect 
to the Maxwell-Boltzmann distribution. Hence, the above equation is well-posed. We 
solve it using separation of variables: 

with 

- £o$ = P- (2.9) 

This Poisson equation is posed on [— 7r,7r] d x M. d . The solution is periodic in q and 
satisfies condition ( 12.61 ). Now we proceed with the third equation in ( 12.81 ). We apply 
the solvability condition to obtain: 



dug 

dt 



— 7T,7T a Xl 



CiUippip, q)dpdq 

d 2 u 



= Y] / Pi$jP0(p,q)dpdq 



This is the Backward Kolmogorov equation which governs the dynamics on large 
scales. We write it in the form 

dun v ^ d 2 Un 

¥ = .^^ (2 ' 10) 

where the effective diffusion tensor is 

Dij = p i $ j p l3 (p,q)dpdq, i,j=l,...d. (2.11) 

J[-ir,Tr] d xM d 

The calculation of the effective diffusion tensor requires the solution of the boundary 
value problem i2.9i and the calculation of the integral in ( 12. 111 ). The limiting backward 
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Kolmogorov equation is well posed since the diffusion tensor is nonnegative. Indeed, 
let (bea unit vector in M. d . We calculate (we use the notation = $ • £ and (-, •) for 
the Euclidean inner product.) 

(£, DO = J(p- W dpdq = J ( - £o* c )$ c /i/9 dpd« 

= IP' 1 J \V P f>i\ 2 ^dpdq>0, (2.12) 

where an integration by parts was used. 

Thus, from the multiscale analysis we conclude that at large lenght/time scales the 
particle which diffuses in a periodic potential performs and effective Brownian motion 
with a nonnegative diffusion tensor which is given by formula ( 12. lit . 

We mention in passing that the analysis presented above can also be applied to 
the problem of Brownian motion in a tilted periodic potential. The Langevin equation 

becomes 

x(t) = -VV(x(t)) +F- 7 i(t) + V '27/3- 1 Wit), (2.13) 

where V(x) is periodic with period 2tt and F is a constant force field. The formulas 
for the effective drift and the effective diffusion tensor are 



V = pp(q,p)dqdp, D = I 

JR d x[-Tr,-ir] d JR d x[-TT,Tr] d 



(p-V)®4>p(p,q)dpdq, (2.14) 



where 



-C<f> = p-V, (2.15a) 
C*p = 0, [ p(p, q) dpdq = 1. (2.15b) 

JR d x\-Tr.ir] d 



with 



C = P ■ V, + (-V,V + F) ■ V p + 7 ( - p ■ V p + /3- x A p ). (2.16) 

We have used ® to denote the tensor product between two vectors; C* denotes the 
L 2 -adjoint of the operator C, i.e. the Fokker-Planck operator. Equations (12.151 are 
equipped with periodic boundary conditions in q. The solution of the Poisson equa- 
tion (12.15b is also taken to be square integrable with respect to the invariant density 

/ l</»(g,p)l 2 p(p, 9) d P d q < +°°- 

JR d x[-n,Tc] d 

The diffusion tensor is nonnegative definite. A calculation similar to the one used to 
derive (I2.12l i shows the positive definiteness of the diffusion tensor: 

(£,D0=lP- 1 J |V P $ ? | 2 pip, q) dpdq >0, (2.17) 

for every vector £ in R d . The study of diffusion in a tilted periodic potential, in the 
underdamped regime and in high dimensions, based on the above formulas for V and 
D, will be the subject of a separate publication. 
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2.2 Equivalence With the Green-Kubo Formula 

Let us now show that the formula for the diffusion tensor obtained in the previous 
section, equation ( 12. lib , is equivalent to the Green-Kubo formula ( 11.41 ). To simplify 
the notation we will prove the equivalence of the two formulas in one dimension. The 
generalization to arbitrary dimensions is immediate. Let (x(t;q,p), v(t;q,p)) with 
v — x and initial conditions x(0;q,p) = q, v(0;q,p) — p be the solution of the 
Langevin equation 

x = —d x V — jx + £ 

where stands for Gaussian white noise in one dimension with correlation function 

(£(*)£(*)) =2 7 k B TS(t-s). 

We assume that the (x, v) process is stationary, i.e. that the initial conditions are 
distributed according to the Maxwell-Boltzmann distribution 

pp{q,p) = Z- 1 e-P H ^\ 

The velocity autocorrelation function is |6l eq. 2.10] 

(v(t;q,p)v(Q;q,p)) = vpp(x,v,t;p,q)pp(p,q)dpdqdxdv, (2.18) 



and p(x, v,t;p, q) is the solution of the Fokker-Planck equation 

= £*P, P(x,v,0;p,q) = 5(x - q)5(v - p), 

where 

C*p = -vd x p + d x V(x)d vP + j{d(vp) + p- l dlp) ■ 
We rewrite ( I2.18l l in the form 

(v(t;q,p)v(0;q,p)) = J J fj J vp(x,v,t;p,q) dvdx^j pppip, q) dpdq 

v{t;p,q)ppp(p,q)dpdq. (2.19) 



The function v(t) satisfies the backward Kolmogorov equation which governs the evo- 
lution of observables I25J Ch. 6] 

-£ = £v, v(0;p,q)=p. (2.20) 
at 

We can write, formally, the solution of ( 12.20b as 

v = e ct p. (2.21) 

We combine now equations ( 12.191 ) and ( 12.211) to obtain the following formula for the 
velocity autocorrelation function 

(v(t;q,p)v(0;q,p)) = [ [ p(e ct p) Pfj (p, q) dpdq. (2.22) 
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We substitute this into the Green-Kubo formula to obtain 



D = 




(v(t;q,p)v{0;q,p))dt 



-L 1 p)pppdpdq 



e dtp I ppp dpdq 



/ / 4>PPi3 dpdq, 



J — OO J — 71 



where <fi is the solution of the Poisson equation (12.91 , In the above derivation we have 
used the formula — Cr 1 = J °° e ct dt, whose proof can be found in 11251 Ch. 11]. 

3 The Underdamped and Overdamped Limits 

In this section we derive approximate formulas for the diffusion coefficient which are 
valid in the overdamped 7^1 and underdampled 7^1 limits. The derivation of 
these formulas is based on the asymptotic analysis of the Poisson equation ( 12.91 ). In 
this section we will take the period of the potential is 2ir. 

3.1 The Underdamped Limit 

In this subsection we solve the Poisson equation j2.9j in one dimension perturbatively 
for small 7. We shall use singular perturbation theory for partial differential equations. 
The operator Cq that appears in (12.9b can be written in the form 



where Ch stands for the (backward) Liouville operator associated with the Hamilto- 
nian H (p, q) and Cqu f° r the generator of the OU process, respectively: 



We expect that the solution of the Poisson equation scales like 7 1 when 7^1. Thus, 
we look for a solution of the form 



Cq = Ch + iCqu 



C H = pd q ~ d q Vd p , Cqu 



-pdp + p^dl 



$ = -00 + 4>\ + 702 + ■ • • 

7 

We substitute this ansatz in j2.9i to obtain the sequence of equations 



(3.1) 



Ch4>q 

Ch4>2 



P + £ou(f>o, 
Cqu4>i- 







(3.2a) 
(3.2b) 
(3.2c) 
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From equation J3.2al > we deduce that, since the 0o is in the null space of the Liouville 
operator, the first term in the expansion is a function of the Hamiltonian z(p,q) = 



\V 2 + V(q): 



0o = <f>o(z(p,q)). 

Now we want to obtain an equation for <f>$ by using the solvability condition for (13.2bl >. 
To this end, we multiply this equation by an arbitrary function of z, g = g(z) and 
integrate over p and q to obtain 



+00 /*7T 



— OO —IT 



(p + C-ouM g(z(p, q)) dpdq = 0. 



We change now fromp, q coordinates to z, q, so that the above integral becomes 



+ OO />7T 



g{z) (p{z,q) + C u(t>o{z)) 



p{z,q) 



dzdq = 0, 



where J — p 1 (z, q) is the Jacobian of the transformation. Operator £0, when applied 
to functions of the Hamiltonian, becomes: 

Cou = (p- 1 -p 2 )-H-+r 1 P 2 -^. 

oz az z 



Hence, the integral equation for (f>o(z) becomes 

+ OO pTT 



g{z) 



p(z, q) + (V 1 ~p 2 )§; + ^Vj^) Mz) 



P(z,q) 



dzdq = 0. 



Let Eq denote the critical energy, i.e. the energy along the separatrix (homoclinic 
orbit). We set 



S(z) 



£Ci(z) 



P(z,q)dq, T(z) 



1 



i(z) 



P(z,q) 



dq, 



where Risken's notation [30j p. 301] has been used for x\{z) and X2(z). 

We need to consider the cases {z > Eq, p > 0} , {z > Eq, p < 0} and {E m i n < 
z < Eq} separately. 

We consider first the case E > Eq, p > 0. In this case x\(x) — tt, X2{z) = —n. 
We can perform the integration with respect to q to obtain 



+00 



9(z) 



2n + UpT^z) - S(z))^- + p^S{z)^j Mz) 



dz = 0, 



This equation is valid for every test function g(z), from which we obtain the following 
differential equation for 0q: 



C<f> := -/?" 



1 

W) 



S(z)<i>" 



1 

W) 



s( z ) - r 1 <p' = 



2tt 
WY 



(3.3) 
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where primes denote differentiation with respect to z and where the subscript has 
been dropped for notational simplicity. 

A similar calculation shows that in the regions E > Eq, p < and E m i n < E < 
Eq the equation for <fio is 

-Zcf> = -^, E>E ,p<0 (3.4) 

and 

-Zcj> = 0, E min <E<E Q . (3.5) 

Equations ( [3.3) , d3.4t . d3.5l > are augmented with condition ( 12.61 > and a continuity con- 
dition at the critical energy (9j 

24,' 3 (E ) = 4>[(E ) + <j>' 2 (E ), (3.6) 

where <px, 02, 4>3 are the solutions of equations ( 13.3b . ( I3.4l i and (13.5b . respectively. 

The average of a function h(q,p) = h(q,p(z, q j) can be written in the form [30] p. 
303] 

/OO />7T 
/ h(q,p)fJ-p{q,p)dqdp 

= W / 

where the partition function is 



(h(q,p(z,q)) + h(q,-p(z,q))j{p(q,z)) 1 e ^ z dzdq, 



Z = x r^-\ e-WMdq. 



P J-* 

From equation (13.51 > we deduce that <j>s{z) = 0. Furthermore, we have that <pi(z) — 
—<p2(z). These facts, together with the above formula for the averaging with respect to 
the Boltzmann distribution, yield: 

D = (p$(p,q))p = (pc(> )p + O(l) (3.7) 
2 r+°° 
« -Zp 1 / M^ z dzO(l) 

7 P ./Eo 

= — Zp 1 / cf> (z)e-^ dz, (3.8) 

to leading order in 7, and where <fio(z) is the solution of the two point boundary value 
problem d3.3l >. We remark that if we start with formula D = 7/3~ 1 (|<9 p $| 2 ) / 3 for the 
diffusion coefficient, we obtain the following formula, which is equivalent to ( 13.81 ): 



\„ r+00 
D=—Zg 1 \d z Mz)\ 2 e~ f3z dz. 
Jeq 



1(3 
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Now we solve the equation for <j>o{z) (f° r notational simplicity, we will drop the sub- 
script ). Using the fact that S'(z) = T(z), we rewrite (13.31 > as 

-/3- 1 (S*0')' + 5*0' = 2tt. 
This equation can be rewritten as 

-0- 1 (e-P*S<l/)=e-f". 
Condition (12. 6t implies that the derivative of the unique solution of (13 - 31 > is 

4>'{z) = s-\z). 



We use this in (13.8b . together with an integration by parts, to obtain the following 
formula for the diffusion coefficient: 

1 f + °° e~$ z 

D = -^ 2Z P 1 ^ 1 "cTTV dz. (3.9) 

7 JEo ^x 2 ) 

We emphasize the fact that this formula is exact in the limit as 7 — > and is valid for 
all periodic potentials and for all values of the temperature. 

Consider now the case of the nonlinear pendulum V(q) — — cos(q). The partition 
function is 

Z - pi/2 J °\Ph 

where Jq(-) is the modified Bessel function of the first kind. Furthermore, a simple 
calculation yields 

S(z) = 2 5/2 Vz~TlE 



z + 1 



where E(-) is the complete elliptic integral of the second kind. The formula for the 
diffusion coefficient becomes 

D = ~ l W T W)li 7fTTEUWTW) dz - (110) 

We use now the asymptotic formula Jq(0) ~ (27r/3) _1 / 2 e' 5 , (3 ^> 1 and the fact that 
E(l) = 1 to obtain the small temperature asymptotics for the diffusion coefficient: 

D=-^-e- 2 ^, /3»1, (3.11) 
7 2/3 

which is precisely formula dl.8t , obtained by Risken. 

Unlike the overdamped limit which is treated in the next section, it is not straight- 
forward to obtain the next order correction in the formula for the effective diffusivity. 
This is because, due to the discontinuity of the solution of the Poisson equation ( 12. 9\ 
along the separatrix. In particular, the next order correction to when 7 -C 1 is of 
(7 -1 / 2 ), rather than (1) as suggested by ansatz (13. It . 
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Upon combining the formula for the diffusion coefficient and the formula for the 
hopping rate from Kramers' theory [ 14, eqn. 4.48(a)] we can obtain a formula for the 
mean square jump length at low friction. For the cosine potential, and for /3 ^> 1, this 
formula is 

<^ 2 > = CT fbr 7 «l,/3»l. (3-12) 



3.2 The Overdamped Limit 

In this subsection we study the large 7 asymptotics of the diffusion coefficient. As in 
the previous case, we use singular perturbation theory, e.g. [ 15 Ch. 8]. The regularity 
of the solution of d2.9t when 7 3> 1 will enable us to obtain the first two terms in the 
- expansion without any difficulty. 

We set 7 = i. The differential operator £ becomes 

£0 = —C-OU + £-H- 

e 

We look for a solution of (12.9b in the form of a power series expansion in 7: 

$ = o + e(/» 1 +e 2 2 +e 3 03 + ... (3.13) 

We substitute this into ( 12. 9\ and obtain the following sequence of equations: 

~C ou ^ = 0, (3.14a) 

-Lou4>\ = p + Ch&o, (3.14b) 

-Coufo = £h(/>i, (3.14c) 

-Coufo = C H (p2- (3.14d) 

The null space of the Ornstein-Uhlenbeck operator Cq consists of constants in p. Con- 
sequently, from the first equation in (13. 14-b we deduce that the first term in the expansion 
in independent of p, <fio = <p(q). The second equation becomes 

-£ou<t>i =p(l + d q 4>). 

Let 

be the invariant distribution of the OU process (i.e. C olJ vp(p) = 0). The solvability 
condition for an equation of the form —Cou4> — f requires that the right hand side 
averages to with respect to v@(p), i.e. that the right hand side of the equation is 
orthogonal to the null space of the adjoint of Cou- This condition is clearly satisfied 
for the equation for <j>i. Thus, by Fredholm alternative, this equation has a solution 
which is 

4>i(p,q) = 0- + dq<t>)p + i>i(q), 
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where the function tpi (q) of is to be determined. We substitute this into the right hand 
side of the third equation to obtain 

-Coufc = P 2 d 2 q - d q V(l + d q cj>) + P d q Mq)- 

From the solvability condition for this we obtain an equation for 4>{q): 

0- 1 d%c/>-d q V(l + d q <t>)=O, (3.15) 

together with the periodic boundary conditions. The derivative of the solution of this 
two-point boundary value problem is 

3 «* +1 =r^? w "- <316) 

The first two terms in the large 7 expansion of the solution of equation d2.9t are 

7 \T 
where 4>(q) is the solution of ( 13.151 l. Substituting this in the formula for the diffusion 
coefficient and using ( 13.16b we obtain 

f°° r 4tt 2 / 1 \ 

d = LJ^M^^ + o^, 

where Z — f n e~^ v ^ q \ Z = e^ v ^ q \ This is, of course, the Lifson-Jackson 
formula which gives the diffusion coefficient in the overdamped limit 1 19 1. Continuing 
in the same fashion, we can also calculate the next two terms in the expansion (13.13b . 
From this, we can compute the next order correction to the diffusion coefficient. The 
final result is 

4tt 2 4ir 2 l3Z 1 m ( 1 



D= '-^ + — , (3.17) 

PjZZ ~f 3 ZZ 2 WJ 

where Z 1 = \V'(q)\ 2 e^ v ^ dq. 

In the case of the nonlinear pendulum, V(q) = cos(q), formula (13.171 1 gives 

D = -W 2 (/3) - 4 f^7§ - Jo 2 (P)) +o(±r), (3.18) 



7/3 " ' 7 3 \W ° ' 7 V7 5 
where J„(/3) is the modified Bessel function of the first kind. 

In the multidimensional case, a similar analysis leads to the large gamma asymp- 
totic s: 

7 \r, 

where £ is an arbitrary unit vector in M. d and Dq is the diffusion coefficient for the 
Smoluchowski (overdamped) dynamics: 



Do = Z- 1 f ( - Cvx) ® xe~ mq) dq 



where 

and x(q) is the solution of the PDE Cvx = V 9 y with periodic boundary conditions. 
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4 The Numerical Method 



In order to calculate the diffusion coefficient D we have to solve equation ( 12.91 ) and 
calculate the integral in ( 12.111 ). It is possible to do this by means of a spectral method. 
In particular, by expanding the solution of ( 12.91 ) into Hermite polynomials in p and a 
standard Fourier series in q, we can convert the boundary value problem into an infinite 
system of linear equations; upon truncating this system we can obtain a finite dimen- 
sional, sparse system of linear equations which we can easily solve. For simplicity we 
will consider the problem in one dimension, though our numerical method works in 
arbitrary dimensions. 

We look for a solution of J2.9I ) in the form 

+oo +oo 
k— — oc n— 

where /„ denotes the nth eigenfunction of the Ornstein-Uhlenbeck operator Cou\ the 
corresponding eigenvalue is k. The eigenfunctions of the OU processes are related to 
the Hermite polynomials through the formula ||30l Ch. 5] 

f n ( P ) = -j=H n (Vj9p) , (4.2) 

where 

«.M -<-ir.*£ («-*)■ 

We remark that the solution of the Poisson equation J2.9I ) is defined up to a constant 
which we have taken to be (i.e. we assume that the solution averages to zero with 
respect to the Maxwell-Boltzmann distribution). It is elementary to check that the 
choice of this constant does not effect the value of the diffusion coefficient. Notice also 
that the boundary conditions (the solution is periodic in q and square integrable with 
respect to the Maxwell-Boltzmann distribution) have already been taken into account 
when writing (14.11 ). 

Using the fact that the p = /3 -1 / 2 /i, together with the orthonormality of the eigen- 
functions of the Ornstein-Uhlenbeck operator, we obtain the following formula for the 
diffusion coefficient: 

OO /-7T 

D = Z^/r 1/2 Yl / e ifc<? e- /3y(9) dq. (4.3) 
For the cosine potential V(q) — — cos(q) the above formula becomes 

D = p- 1,2 Ja\P) fakMP), (4.4) 

k— — oo 

Now we shall obtain a linear system of equations for the coefficients <fi n k for the cosine 
potential. The Poisson equation ( |2.9t in one dimension is 

-£<P = p = f3- 1/2 A- (4.5) 
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1(T 3 1(T 2 10"' 10° 10 1 10 2 10 3 

Y 

Figure 1 : Diffusion coefficient as a function of the friction coefficient for the cosine 
potential. Dash-dot line and circles: D obtained from the numerical solution of the 
Poisson equation, formula (|4.4| i; stars: D obtained from the the calculation of the jump 
length distribution and the hopping rate, formula dl.2| ):diamonds: results from Monte 
Carlo simulations, formula ( 1 1 .3b ;solid lines, analytical approximation for 7 -C 1, 7 ^ 
1, equations ( 13. lit and ( 11.7b . 



We substitute the expansion ( 14. It into ( 14. 5\ and use the properties of Hermite polyno- 
mials and of trigonometric functions to obtain the following linear system of equations 

2[3- 1 kc/) lk + cl) lk ^ 1 -<j> lk+1 =0, (4.6a) 

2 7 z0! k + 2 y/^ktpa ft + 2 v / 2/F T fc0 2 fe + ^2/302 fe-i - y/20<h k+i = ^ ^JW 1 4 , 

(4.6b) 

1\J fi~ x \Jn + lkcj) n+ i k + 2\f^-^^/nk4>n-ik 

+ y^Vn+ l<f> n+lk -i - v^V n + l^n+i fc+i + 2i-fntf> nk = 0. (4.6c) 

We truncate ( I4.6I > by taking into account the first N + 1 terms of the Hermite expansion 
and the first K + 1 terms of the Fourier expansion. This leads to a sparse linear system 
of (N + 1) (2K + 1) equations and (N + 1) (2K + 1) - 1 variables (note the absence 
of the term <fioo). This implies that one of the equations in ( 14. 6t is linearly dependent 
on the others; to obtain a nonsingular system we remove the equation for k — 00 In 

3 This follows from our assumption that the solution of the Poisson equation averages to with respect to 
the canonical ensemble. 
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this way we obtain a nonsingular system of s = (N + 1) (2K + 1) — 1 equations and 
unknowns. This system can be written in the form Ax = b with 



$n k 7 i 



k + K + 1 n = 0,k<0, 

n (2K + 1) + k + K n = 0, k > or n ^ 



and otij is the coefficient of Xj in the i th equation, which is taken from equations ( 14. 6t 
for values of n and k 



The sparsity of this system implies that we can solve it very efficiently. In particular, 
we can calculate accurately the diffusion coefficient with a minimal computational cost, 
even for very small values of the friction constant 7. 

To illustrate the efficiency of our numerical method, we calculate the diffusion co- 
efficient for the cosine potential as a function of the dissipation 7, at a fixed temperature 
/3 _1 = 0.5. We compare the results obtained through our numerical method with the 
approximate analytical expressions dl.7| > and dl.81 l, results from Monte Cairo simu- 
lations using d 1 . 31 > and results obtained through numerical calculation of the hopping 
rate and the jump length distribution, equation ( 11.21 ). FOr the calculation of (£ 2 ) and 
k we generate a long path (for every 7) of the Langevin dynamics using the Milstein 
scheme. 

The results of the numerical simulations are presented in Figure [4] There agree- 
ment between the approximate analytical formulas and the calculation of the diffusion 
coefficient using the method described in this section are excellent. Our method is 
far superior in comparison to Monte Carlo simulations or the calculation of the mean 
square jump length and the hopping rate, since even for very small 7 the solution of a 
rather small linear system of equations is required^ On the other hand, contrast, the 
path of integration over which we have to solve the Langevin equation in order to com- 
pute accurate statistics increases as 7 decreases and the calculation of D using Monte 
Carlo becomes computationally expensive. 

5 Conclusions 

The problem of Brownian motion in a periodic potential in arbitrary dimensions was 
studied in this paper. Using multiscale techniques ll25l we derived a formula for the 
effective diffusion tensor which is valid for all values of the friction coefficient and the 
temperature, and in arbitrary dimensions. We also showed the equivalence between 
our formula and the Green-Kubo formula. The calculation of the diffusion tensor using 
our approach requires the solution of a Poisson equation together with the calculation 
of the average of an appropriate function with respect to the canonical distribution. 
Furthermore, the overdamped and underdamped limits where studied and approximate 

4 Needless to say, more Hermite and Fourier terms have to be taken into account when 7 decreases. 
However, even for 7 very small, the resulting lineal' system of equations is small enough so that it can be 
solved in a few seconds in Matlab. 
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analytical formulas for the diffusion coefficient in these two limits where obtained. In 
addition, a very efficient numerical method for the calculation of the diffusion coef- 
ficient was developed; this numerical method is based on the solution of the Poisson 
equation via a spectral method and it leads to the accurate and very efficient calculation 
of the diffusion coefficient even for very low values of the friction coefficient. 

The approach developed in the paper for the study of the problem of diffusion 
in periodic potentials offers various advantages over other analytical and numerical 
methods. First, all the results reported in this paper can be justified rigorously and 
they can also lead to a rigorous analysis of the dependence of the diffusion tensor 
on the friction coefficient and on the temperature. A first step in this direction was 
taken in JT2). Second, our method enables us to study various distinguished limits 
of physical interest (such as the overdamped and underdamped limits) in a systematic 
fashion through asymptotic analysis of the Poisson equation. Third, it leads to an 
efficient numerical method for calculating the diffusion tensor through the numerical 
solution of the Poisson equation. The effectiveness of our method was shown in this 
paper for the one dimensional problem, for which analytical approximate formulas are 
can be derived. The method is also very efficient in two and three dimensions and 
can offer insight into the problem of Brownian motion in periodic potentials in higher 
dimensions. A thorough numerical investigation of the multidimensional problem will 
be presented elsewhere. 
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